home *** CD-ROM | disk | FTP | other *** search
/ FM Towns: Free Software Collection 4 / FM Towns Free Software Collection 4 - Disc 1.iso / t_os / pi_v2 / pi.c next >
Encoding:
C/C++ Source or Header  |  1991-10-19  |  32.2 KB  |  946 lines
  1. /* ANSI-C(?)によるπ計算   VERSION 2.0 */
  2. /* Copyright(c) Daisuke Takahashi 1991 */
  3.  
  4. /* 任意の瞬間にキー入力によって計算を中断したいのですが */
  5. /* High-Cで signal.h が使えないので、一定時間毎にセーブさせています */
  6. /* その代わり、DOS.H未使用の為、殆どのＣでコンパイルできます(^_^) */
  7.  
  8. #include <stdio.h>
  9. #include <stdlib.h>
  10. #include <ctype.h>
  11. #include <time.h>
  12. #include <string.h>
  13.  
  14. /* ----------- definitions ------------ */
  15. #define TEN8   100000000L  /* １億 10^8 */
  16. #define MAX    4190000000L /* max of USL - TEN8 */
  17. #define ON     (1)
  18. #define OFF    (0)
  19. #define SET    (2)
  20. #define LAP    (1)
  21. #define DONE   (0)
  22. #define DATA   (1)
  23. #define PLUS   (1)
  24. #define MINUS  (-1)
  25. #define NUL    (\0)
  26. #define USL    unsigned long     /* typedef にすべきでしょうが(^_^;) */
  27. #define RATE   (int)( 1000 * (long)ex.head / (long)ex.max ) 
  28.  
  29. /* ---------------- prototypes ------------------ */
  30.   void set_ex(int,char**);              /* 引数→変数 */
  31.   void get_memory(void);
  32.   void back_memory(void);               /* ワーク返還 */
  33.   void formal(void);                    /* 計算公式   */
  34.   void arctan(int,int,int);             /* arctan設定 */
  35.   long bench_mark(long);                /* 速度測定   */
  36.   void message(int,int,int);            /* 計算量表示 */
  37.   void c_std(int,USL,USL,long*,long*);  /* 41以下     */
  38.   void c_low(int,USL,USL,long*,long*);  /* 6.4万以下  */
  39.   void c_mid(int,USL,USL,long*,long*);  /* 200万以下  */
  40.   void c_high(int,USL,USL,long*,long*); /* 256万以下  */
  41.   void c_huge(int,USL,USL,long*,long*); /* 1600万以下 */
  42.   void regular(int);                    /* 正規化     */
  43.   void load(char *);                    /* 経過ロード */
  44.   void save(void);                      /* 経過セーブ */
  45.   void show_ans(void);                  /* 結果の表示 */
  46.   void time_set(void);                  /* 開始時刻記録 */
  47.   void time_end(void);                  /* 計算時間更新 */
  48.   void usage(int);                      /* 使い方 */
  49.  
  50. /* ------------- structs --------- */
  51. struct {
  52.   /* ファイルセーブする変数 */
  53.   int  ver;   /* Version  */
  54.   long sec;   /* 計算時間 */
  55.   long cent;  /* 1/100秒  */
  56.   long keta;  /* 計算桁数 */
  57.   int  max;   /* 配列数   */
  58.   int  type;  /* 計算公式 */
  59.   int  head;  /* 動的先頭 */
  60.   int  tail;  /* 動的最終 */
  61.   USL  save;  /* 退避flag */
  62.   USL  con;   /*          */
  63.   USL  var;   /*          */
  64.   /* セーブ必要なし */
  65.   long *arc;      /* 計算作業領域 */
  66.   long *ans;      /* 計算結果領域 */
  67.   int  load;      /* file読込済   */
  68.   time_t  s_sec;  /* 計算開始時刻 */
  69.   clock_t s_cent; /* 〃1/100 秒   */
  70. } ex;
  71.  
  72. int main( int argc, char **argv ) {
  73.   set_ex( argc, argv );  /* 引き数から変数設定 */
  74.  
  75.   formal();      /* 計算 */
  76.  
  77.   show_ans();    /* 結果と時間を表示 */
  78.   back_memory(); /* ワーク返還 */
  79.   return 0;      /* 正常終了   */
  80. }
  81.  
  82. void set_ex( int argc, char **argv ) {
  83.   char *buf;
  84.  
  85.   if (argc == 1)                       /* オプション無しなら */
  86.     usage( !DATA );                  /* 普通説明表示       */
  87.   if (strchr( argv[1], '.' ) != NULL)  /* ファイル名であれば */
  88.     load( argv[1] );                 /* 途中結果を読み込む */
  89.   else if ( !(isdigit( *argv[1] ) ) )  /* 無効オプションなら */
  90.     usage( DATA );                   /* データ説明         */
  91.   else {            /* 桁数指定オプション  */
  92.     ex.ver  = 0;  /* On memory Version   */
  93.     ex.sec  = 0;  /* 計算時間(秒単位)    */
  94.     ex.cent = 0;  /* 計算時間(0.01秒単位 */
  95.     ex.keta = atol( argv[1] );         /* 計算桁を代入 */
  96.     ex.max  = (int)(ex.keta / 8 + 2);  /* 配列数 */
  97.     ex.save = 0;  /* セーブモードでない */
  98.     get_memory(); /* 必要メモリを確保   */
  99.   }
  100.   if (argc == 3) {
  101.     buf = argv[2];
  102.     while ( *buf ) {
  103.       if ( isdigit(*buf) )                    /* 数字の場合は */
  104.         ex.save = strtol( buf, &buf, 10 );  /* 分単位のセーブ時間 */
  105.       else {
  106.         if (strchr( "gmrsGMRS", *buf ) != NULL && ex.load == OFF)
  107.           ex.type = *buf;  /* ロード後の公式変更は不可 */
  108.         buf++;               /* 計算公式を代入 */
  109.       }
  110.     }
  111.   }
  112. }
  113.  
  114. void get_memory( void ) {      /* ワークを確保 */
  115.   ex.arc = (long *)calloc( (size_t)ex.max, (size_t)sizeof(long) ); 
  116.   ex.ans = (long *)calloc( (size_t)ex.max, (size_t)sizeof(long) );
  117.   if (ex.ans == NULL) {
  118.     fputs( "メモリーが足らんぜよ(;_;) \n", stderr );
  119.     exit( 1 );
  120.   }
  121. }
  122.  
  123. void back_memory( void ) {   /* ワークを返還 */
  124.   free( ex.arc );
  125.   free( ex.ans );
  126. }
  127.  
  128. void formal( void ) {    /* 計算公式の値をセット */
  129.   switch ( tolower( ex.type ) ) {
  130.  
  131.     case 'g': /* ガウスの式 */
  132.       /* 48*arctan(1/18) + 32*arctan(1/57) - 20*arctan(1/239) */
  133.       arctan( PLUS, 48, 18 );
  134.       arctan( PLUS, 32, 57 );
  135.       arctan( MINUS, 20, 239 );
  136.       break;
  137.     case 'r': /* ラザフォードの式 */
  138.       /* 16*arctan(1/5) - 4*arctan(1/70) + 4*arctan(1/99) */
  139.       arctan( PLUS, 16, 5 );
  140.       arctan( MINUS, 4, 70 );
  141.       arctan( PLUS, 4, 99 );
  142.       break;
  143.     case 's': /* シュテルマーの式 */
  144.       /* 24*arctan(1/8) + 8*arctan(1/57) + 4*arctan(1/239) */
  145.       arctan( PLUS, 24, 8 );
  146.       arctan( PLUS, 8, 57 );
  147.       arctan( PLUS, 4, 239 );
  148.       break;
  149.     default:  /* デフォルトは最速マチン式 */
  150.     case 'm': /* マチンの式 */
  151.       /* 16*arctan(1/5) - 4*arctan(1/239) */
  152.       arctan( PLUS, 16, 5 );
  153.       arctan( MINUS, 4, 239 );
  154.       break;
  155.   }
  156.   if (ex.load == ON) 
  157.     fputs( "Not calculated !!!\n", stderr );
  158. }
  159.  
  160. void arctan( int sign, int mag, int base ) {
  161.   /* arctan(1/base) = 1/base - 1/(3*base^3) + 1/(5*base^5) - */
  162.   time_t cur_time, pre_time;
  163.   long min_one;     /* １分相当の var */
  164.  
  165.   message( SET, sign, SET );  /* 3番目はダミー：計算式の符号 */
  166.   if (ex.load == OFF ) {    /* 新規計算開始 */
  167.     *ex.arc = (long)base * (long)mag;  /* 初期値設定 */
  168.     ex.con  = (long)base * (long)base;
  169.     ex.head = 0;
  170.     ex.var  = 1;
  171.     if ( base == 5 || base == 8 )  /* 10進数で有限小数なら */
  172.       ex.tail = 1;                 /* 最終位置は動的に変化 */
  173.     else                 /* 10進数で循環小数なら最終位置は物理的最終 */
  174.       ex.tail = ex.max;
  175.   } else {                  /* 計算再開 */
  176.     if (ex.con != (long)base * (long)base) {
  177.       message( DONE, mag, base );
  178.       return;  /* この項は計算済 */
  179.     } else 
  180.       ex.load = OFF;   /* この項から計算再開：ロードフラグクリア */
  181.   }
  182.   message( LAP, mag, base );
  183.  
  184.   time_set();   /* ここから計算時間測定開始 */
  185.   time( &pre_time );  /* 表示間隔初期値 */
  186.   min_one = bench_mark( 0 );  /* 時間測定初期化 */
  187.  
  188.   while ( (ex.head < ex.max) || *(ex.arc + ex.max) ) {
  189.     USL    check;       /* どのcalc_*_? を使うかの判定 */
  190.  
  191.     if ( ex.var > 64000L && (TEN8 % ex.var) < (MAX / ex.var) ) { 
  192.       check = 60000L;     /* enough under 64000L */
  193.       /* ans_r * T8_v_r < 41億 なら calc_low() で計算 */
  194.     } else 
  195.       check = ex.var;
  196.     if      (check < 42)
  197.       c_std( sign, ex.con, ex.var, ex.arc, ex.ans ); 
  198.     else if (check < 64000L) 
  199.       c_low( sign, ex.con, ex.var, ex.arc, ex.ans );
  200.     else if (check < 2030000L)
  201.       c_mid( sign, ex.con, ex.var, ex.arc, ex.ans );
  202.     else if (check < 2560000L)
  203.       c_high( sign, ex.con, ex.var, ex.arc, ex.ans );
  204.     else /*  check < 1600万 */
  205.       c_huge( sign, ex.con, ex.var, ex.arc, ex.ans );
  206.  
  207.     sign *= MINUS; /* 符号逆転 */
  208.       if ( ex.tail < ex.max && *(ex.arc + ex.tail) )
  209.         ex.tail++;
  210.       if ( !*(ex.arc + ex.head) && ex.head < ex.max ) 
  211.         ex.head++;  /* 有効数字の先頭に移動 */
  212.       ex.var += 2;  /* 次の展開項を計算     */
  213.  
  214.     if ( ex.var % 80 == 1 ) { /* ex.var 80  毎に */
  215.       regular( OFF );       /* ans[]を正規表現化（高速モード） */
  216.  
  217.       if (min_one < 80)    /* １分相当値未計算なら */
  218.         min_one = bench_mark( min_one );  /* 計算する */
  219.  
  220.       if ( (min_one >= 80) && ex.var % min_one == 1 ) {  /* １分経過なら */
  221.         time( &cur_time );  /* 現在時刻の取得 */
  222.  
  223.         if (ex.save) {   /* セーブ指定有りで時間なら */
  224.           if ( (long)difftime( cur_time, ex.s_sec ) > 60 * ex.save )
  225.             save();      /* セーブするだよん♪(^0^) */
  226.         }
  227.         if ( (long)difftime( cur_time, pre_time ) <= 40 ) {
  228.           min_one *= 2;   /* 表示間隔が短くなったら倍にする */
  229.         }
  230.         pre_time = cur_time;
  231.         message( LAP, mag, base ); /* 表示時 Ctrl-C有効 */
  232.       }
  233.     }
  234.   }
  235.   regular( ON );
  236.   time_end();    /* 次のarctan()の為 */
  237.   message( DONE, mag, base );
  238.   if (ex.save)  /* 一つのarctan(1/n)終了時も */
  239.     save();     /* saveする */
  240. }
  241.  
  242. long bench_mark( long min_one ) {   /* 1 分相当の ex.var 差を返す */
  243.   time_t cur_time;
  244.   long   sec;
  245.   static USL st_var;
  246.  
  247.   if (min_one == 0) { /* 初回計測？   */
  248.     st_var = ex.var;  /* 計測開始時刻 */
  249.     return 1;         /* 初回時は計測 */
  250.   }
  251.   if (min_one >= 2)   /* まだまだ時間がある */
  252.     return min_one - 1;
  253.  
  254.   time( &cur_time );
  255.   sec = (long)difftime( cur_time, ex.s_sec );
  256.   if (sec >= 60)
  257.     return  ex.var - st_var; /* １分相当値：多桁のロスを防ぐ為80加算 */
  258.   else if (sec >= 30)       /* 80加算の理由は上と同じ */
  259.     return (ex.var - st_var) *60 /sec /80 *80;  /* 1分近似値 */
  260.   else                      /* 30秒以上余裕がある場合は */
  261.     return (60 - sec) / (sec + 1);  /* 残り回数の推算 */
  262. }
  263.  
  264. void message( int mode, int mag, int base ) {
  265.   static int sign;
  266.  
  267.   if (mode == SET)
  268.     sign = mag;  /* 計算式全体のフラグ */
  269.  
  270.   if (mode == LAP)
  271.     fprintf( stderr, "%darctan(1/%d) %2d.%d%%\r",
  272.                       sign * mag, base, RATE/10, RATE%10 );
  273.   if (mode == DONE)
  274.     fprintf( stderr, "%darctan(1/%d) completed !!!\n", sign * mag, base );
  275. }
  276.  
  277. void c_std( int sign, USL con, USL var, long *arc, long *ans ) {
  278.   /* var must be under 42 */
  279.   int i, j;
  280.   USL co, va, buf, arc_r, ans_r;
  281.   USL con_q = TEN8 / con;  /* arctan(1/5)以外 */
  282.   USL con_r = TEN8 % con;  /* arctan(1/5)以外 */
  283.  
  284.   if (sign == PLUS) {
  285.     if (con <= 40) {   /* arctan(1/5)のみ */
  286.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  287.            i = ex.head, j = ex.tail; i <= j; i++) {
  288.         buf      = *(arc+i) + arc_r * TEN8;
  289.         *(arc+i) = buf / co;
  290.         arc_r    = buf % co;
  291.         buf       = *(arc+i) + ans_r * TEN8;
  292.         *(ans+i) += buf / va;
  293.         ans_r     = buf % va;
  294.       }
  295.       for (j = ex.max; i <= j; i++) {
  296.         buf       = ans_r * TEN8;
  297.         *(ans+i) += buf / va;
  298.         ans_r     = buf % va;
  299.       }
  300.     } else {         /* arctan(1/5)以外 */
  301.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  302.            i = ex.head, j = ex.tail; i <= j; i++) {
  303.         buf      = *(arc+i) + arc_r * con_r;
  304.         *(arc+i) = buf / co + arc_r * con_q;
  305.         arc_r    = buf % co;
  306.         buf       = *(arc+i) + ans_r * TEN8;
  307.         *(ans+i) += buf / va;
  308.         ans_r     = buf % va;
  309.       }
  310.       for (j = ex.max; i <= j; i++) {
  311.         buf       = ans_r * TEN8;
  312.         *(ans+i) += buf / va;
  313.         ans_r     = buf % va;
  314.       }
  315.     }
  316.   } else {    /* sign == MINUS */
  317.     if (con <= 40) {  /* arctan(1/5)のみ */
  318.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  319.            i = ex.head, j = ex.tail; i <= j; i++) {
  320.         buf      = *(arc+i) + arc_r * TEN8;
  321.         *(arc+i) = buf / co;
  322.         arc_r    = buf % co;
  323.         buf       = *(arc+i) + ans_r * TEN8;
  324.         *(ans+i) -= buf / va;
  325.         ans_r     = buf % va;
  326.       }
  327.       for (j = ex.max; i <= j; i++) {
  328.         buf       = ans_r * TEN8;
  329.         *(ans+i) -= buf / va;
  330.         ans_r     = buf % va;
  331.       }
  332.     } else {         /* arctan(1/5)以外 */
  333.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  334.            i = ex.head, j = ex.tail; i <= j; i++) {
  335.         buf      = *(arc+i) + arc_r * con_r;
  336.         *(arc+i) = buf / co + arc_r * con_q;
  337.         arc_r    = buf % co;
  338.         buf       = *(arc+i) + ans_r * TEN8;
  339.         *(ans+i) -= buf / va;
  340.         ans_r     = buf % va;
  341.       }
  342.       for (j = ex.max; i <= j; i++) {
  343.         buf       = ans_r * TEN8;
  344.         *(ans+i) -= buf / va;
  345.         ans_r     = buf % va;
  346.       }
  347.     }
  348.   }
  349. }
  350.  
  351. void c_low( int sign, USL con, USL var, long *arc, long *ans ) {
  352.   /* var must be under 64_000 */
  353.   int i, j;
  354.   USL co, va, buf, arc_r, ans_r;
  355.   USL con_q = TEN8 / con;  /* arctan(1/5)以外 */
  356.   USL con_r = TEN8 % con;  /* arctan(1/5)以外 */
  357.   USL var_q = TEN8 / var;
  358.   USL var_r = TEN8 % var;
  359.  
  360.   if (sign == PLUS) {
  361.     if (con <= 40) {  /* arctan(1/5)のみ */
  362.       for ( co = con, va = var, arc_r = 0, ans_r = 0,
  363.             i = ex.head, j = ex.tail; i <= j; i++) {
  364.         buf      = *(arc+i) + arc_r * TEN8;
  365.         *(arc+i) = buf / co;
  366.         arc_r    = buf % co;
  367.         buf       = *(arc+i) + ans_r * var_r;
  368.         *(ans+i) += buf / va + ans_r * var_q;
  369.         ans_r     = buf % va;
  370.       }
  371.       for (j = ex.max; i <= j; i++ ) {
  372.         buf       = ans_r * var_r;
  373.         *(ans+i) += buf / va + ans_r * var_q;
  374.         ans_r     = buf % va;
  375.       }
  376.     } else {         /* arctan(1/5)以外 */
  377.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  378.            i = ex.head, j = ex.tail; i <= j; i++) {
  379.         buf      = *(arc+i) + arc_r * con_r;
  380.         *(arc+i) = buf / co + arc_r * con_q;
  381.         arc_r    = buf % co;
  382.         buf       = *(arc+i) + ans_r * var_r;
  383.         *(ans+i) += buf / va + ans_r * var_q;
  384.         ans_r     = buf % va;
  385.       }
  386.       for (j = ex.max; i <= j; i++ ) {
  387.         buf       = ans_r * var_r;
  388.         *(ans+i) += buf / va + ans_r * var_q;
  389.         ans_r     = buf % va;
  390.       }
  391.     }
  392.   } else {    /* sign == MINUS */
  393.     if (con <= 40) {  /* arctan(1/5)のみ */
  394.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  395.            i = ex.head, j = ex.tail; i <= j; i++) {
  396.         buf      = *(arc+i) + arc_r * TEN8;
  397.         *(arc+i) = buf / co;
  398.         arc_r    = buf % co;
  399.         buf       = *(arc+i) + ans_r * var_r;
  400.         *(ans+i) -= buf / va + ans_r * var_q;
  401.         ans_r     = buf % va;
  402.       }
  403.       for (j = ex.max; i <= j; i++ ) {
  404.         buf       = ans_r * var_r;
  405.         *(ans+i) -= buf / va + ans_r * var_q;
  406.         ans_r     = buf % va;
  407.       }
  408.     } else {         /* arctan(1/5)以外 */
  409.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  410.            i = ex.head, j = ex.tail; i <= j; i++) {
  411.         buf      = *(arc+i) + arc_r * con_r;
  412.         *(arc+i) = buf / co + arc_r * con_q;
  413.         arc_r    = buf % co;
  414.         buf       = *(arc+i) + ans_r * var_r;
  415.         *(ans+i) -= buf / va + ans_r * var_q;
  416.         ans_r     = buf % va;
  417.       }
  418.       for (j = ex.max; i <= j; i++ ) {
  419.         buf       = ans_r * var_r;
  420.         *(ans+i) -= buf / va + ans_r * var_q;
  421.         ans_r     = buf % va;
  422.       }
  423.     }
  424.   }
  425. }
  426.  
  427. void c_mid( int sign, USL con, USL var, long *arc, long *ans ) {
  428.   /* var must be under 2_030_000 */
  429.   int i, j;
  430.   USL co, va, buf, buf2, quot;
  431.   USL arc_r, ans_r;
  432.   USL con_q  = TEN8 / con;  /* arctan(1/5)以外 */
  433.   USL con_r  = TEN8 % con;  /* arctan(1/5)以外 */
  434.   USL var_q  = TEN8 / var;
  435.   USL var_r1 = (TEN8 % var) / 1024;
  436.   USL var_r2 = (TEN8 % var) % 1024;
  437.  
  438.   if (sign == PLUS) {
  439.     if (con <= 40) {  /* arctan(1/5)のみ */
  440.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  441.            i = ex.head, j = ex.tail; i <= j; i++) {
  442.         buf      = *(arc+i) + arc_r * TEN8;
  443.         arc_r    = buf % co;
  444.         *(arc+i) = buf / co;
  445.         buf    = ans_r * var_r1;
  446.         quot   = ( (buf / va) << 10 ) + ans_r * var_q;
  447.         buf2   = ( (buf % va) << 10 ) + ans_r * var_r2 + *(arc+i);
  448.         *(ans+i) += buf2 / va + quot;
  449.         ans_r     = buf2 % va;
  450.       }
  451.       for (j = ex.max; i <= j; i++) {
  452.         buf    = ans_r * var_r1;
  453.         quot   = ( (buf / va) << 10 ) + ans_r * var_q;
  454.         buf2   = ( (buf % va) << 10 ) + ans_r * var_r2;
  455.         *(ans+i) += buf2 / va + quot;
  456.         ans_r     = buf2 % va;
  457.       }
  458.     } else {         /* arctan(1/5)以外 */
  459.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  460.            i = ex.head, j = ex.tail; i <= j; i++) {
  461.         buf      = *(arc+i) + arc_r * con_r;
  462.         *(arc+i) = buf / co + arc_r * con_q;
  463.         arc_r    = buf % co;
  464.         buf    = ans_r * var_r1;
  465.         quot   = ( (buf / va) << 10 ) + ans_r * var_q;
  466.         buf2   = ( (buf % va) << 10 ) + ans_r * var_r2 + *(arc+i);
  467.         *(ans+i) += buf2 / va + quot;
  468.         ans_r     = buf2 % va;
  469.       }
  470.       for (j = ex.max; i <= j; i++) {
  471.         buf    = ans_r * var_r1;
  472.         quot   = ( (buf / va) << 10 ) + ans_r * var_q;
  473.         buf2   = ( (buf % va) << 10 ) + ans_r * var_r2;
  474.         *(ans+i) += buf2 / va + quot;
  475.         ans_r     = buf2 % va;
  476.       }
  477.     }
  478.   } else {    /* sign == MINUS */
  479.     if (con <= 40) {  /* arctan(1/5)のみ */
  480.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  481.            i = ex.head, j = ex.tail; i <= j; i++) {
  482.         buf      = *(arc+i) + arc_r * TEN8;
  483.         arc_r    = buf % co;
  484.         *(arc+i) = buf / co;
  485.         buf    = ans_r * var_r1;
  486.         quot   = ( (buf / va) << 10 ) + ans_r * var_q;
  487.         buf2   = ( (buf % va) << 10 ) + ans_r * var_r2 + *(arc+i);
  488.         *(ans+i) -= buf2 / va + quot;
  489.         ans_r     = buf2 % va;
  490.       }
  491.       for (j = ex.max; i <= j; i++) {
  492.         buf    = ans_r * var_r1;
  493.         quot   = ( (buf / va) << 10 ) + ans_r * var_q;
  494.         buf2   = ( (buf % va) << 10 ) + ans_r * var_r2;
  495.         *(ans+i) -= buf2 / va + quot;
  496.         ans_r     = buf2 % va;
  497.       }
  498.     } else {         /* arctan(1/5)以外 */
  499.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  500.            i = ex.head, j = ex.tail; i <= j; i++) {
  501.         buf      = *(arc+i) + arc_r * con_r;
  502.         *(arc+i) = buf / co + arc_r * con_q;
  503.         arc_r    = buf % co;
  504.         buf    = ans_r * var_r1;
  505.         quot   = ( (buf / va) << 10 ) + ans_r * var_q;
  506.         buf2   = ( (buf % va) << 10 ) + ans_r * var_r2 + *(arc+i);
  507.         *(ans+i) -= buf2 / va + quot;
  508.         ans_r     = buf2 % va;
  509.       }
  510.       for (j = ex.max; i <= j; i++) {
  511.         buf    = ans_r * var_r1;
  512.         quot   = ( (buf / va) << 10 ) + ans_r * var_q;
  513.         buf2   = ( (buf % va) << 10 ) + ans_r * var_r2;
  514.         *(ans+i) -= buf2 / va + quot;
  515.         ans_r     = buf2 % va;
  516.       }
  517.     }
  518.   }
  519. }
  520.  
  521. void c_high( int sign, USL con, USL var, long *arc, long *ans ) {
  522.   /* var must be under 2_560_000 */
  523.   int i, j;
  524.   USL co, va, buf, buf2, quot;
  525.   USL arc_r, ans_r;
  526.   USL con_q  = TEN8 / con;  /* arctan(1/5)以外 */
  527.   USL con_r  = TEN8 % con;  /* arctan(1/5)以外 */
  528.   USL var_q  = TEN8 / var;
  529.   USL var_r1 = (TEN8 % var) / 1600;
  530.   USL var_r2 = (TEN8 % var) % 1600;
  531.  
  532.   if (sign == PLUS) {
  533.     if (con <= 40) {
  534.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  535.            i = ex.head, j = ex.tail; i <= j; i++) {
  536.         buf      = *(arc+i) + arc_r * TEN8;
  537.         *(arc+i) = buf / co;
  538.         arc_r    = buf % co;
  539.         buf    = ans_r * var_r1;
  540.         quot   = ans_r * var_q;
  541.         quot  += 1600 * (buf / va);
  542.         buf2   = 1600 * (buf % va) + ans_r * var_r2 + *(arc+i);
  543.         *(ans+i) += buf2 / va + quot;
  544.         ans_r     = buf2 % va;
  545.       }
  546.       for (j = ex.max; i <= j; i++) {
  547.         buf    = ans_r * var_r1;
  548.         quot   = ans_r * var_q;
  549.         quot  += 1600 * (buf / va);
  550.         buf2   = 1600 * (buf % va) + ans_r * var_r2;
  551.         *(ans+i) += buf2 / va + quot;
  552.         ans_r     = buf2 % va;
  553.       }
  554.     } else {       
  555.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  556.            i = ex.head, j = ex.tail; i <= j; i++) {
  557.         buf      = *(arc+i) + arc_r * con_r;
  558.         *(arc+i) = buf / co + arc_r * con_q;
  559.         arc_r    = buf % co;
  560.         buf    = ans_r * var_r1;
  561.         quot   = ans_r * var_q;
  562.         quot  += 1600*(buf / va);
  563.         buf2   = 1600*(buf % va) + ans_r * var_r2 + *(arc+i);
  564.         *(ans+i) += buf2 / va + quot;
  565.         ans_r     = buf2 % va;
  566.       }
  567.       for (j = ex.max; i <= j; i++) {
  568.         buf    = ans_r * var_r1;
  569.         quot   = ans_r * var_q;
  570.         quot  += 1600 * (buf / va);
  571.         buf2   = 1600 * (buf % va) + ans_r * var_r2;
  572.         *(ans+i) += buf2 / va + quot;
  573.         ans_r     = buf2 % va;
  574.       }
  575.     }
  576.   } else {
  577.     if (con <= 40) {
  578.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  579.            i = ex.head, j = ex.tail; i <= j; i++) {
  580.         buf      = *(arc+i) + arc_r * TEN8;
  581.         *(arc+i) = buf / co;
  582.         arc_r    = buf % co;
  583.         buf    = ans_r * var_r1;
  584.         quot   = ans_r * var_q;
  585.         quot  += 1600 * (buf / va);
  586.         buf2   = 1600 * (buf % va) + ans_r * var_r2 + *(arc+i);
  587.         *(ans+i) -= buf2 / va + quot;
  588.         ans_r     = buf2 % va;
  589.       }
  590.       for (j = ex.max; i <= j; i++) {
  591.         buf    = ans_r * var_r1;
  592.         quot   = ans_r * var_q;
  593.         quot  += 1600 * (buf / va);
  594.         buf2   = 1600 * (buf % va) + ans_r * var_r2;
  595.         *(ans+i) -= buf2 / va + quot;
  596.         ans_r     = buf2 % va;
  597.       }
  598.     } else {
  599.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  600.            i = ex.head, j = ex.tail; i <= j; i++) {
  601.         buf      = *(arc+i) + arc_r * con_r;
  602.         *(arc+i) = buf / co + arc_r * con_q;
  603.         arc_r    = buf % co;
  604.         buf    = ans_r * var_r1;
  605.         quot   = ans_r * var_q;
  606.         quot  += 1600 * (buf / va);
  607.         buf2   = 1600 * (buf % va) + ans_r * var_r2 + *(arc+ i);
  608.         *(ans+i) -= buf2 / va + quot;
  609.         ans_r     = buf2 % va;
  610.       }
  611.       for (j = ex.max; i <= j; i++) {
  612.         buf    = ans_r * var_r1;
  613.         quot   = ans_r * var_q;
  614.         quot  += 1600 * (buf / va);
  615.         buf2   = 1600 * (buf % va) + ans_r * var_r2;
  616.         *(ans+i) -= buf2 / va + quot;
  617.         ans_r     = buf2 % va;
  618.       }
  619.     }
  620.   }
  621. }
  622.  
  623. void c_huge( int sign, USL con, USL var, long *arc, long *ans ) {
  624.   /* var must be under 15_625_000 (=250^3) */
  625.   int i, j;
  626.   USL co, va, buf, buf2, quot;
  627.   USL arc_r, ans_r;
  628.   USL con_q = TEN8 / con;  /* arctan(1/5)以外 */
  629.   USL con_r = TEN8 % con;  /* arctan(1/5)以外 */
  630.   USL var_q = TEN8 / var;
  631.   USL var_r1 = (TEN8 % var) / 62500L;
  632.   USL var_r2 = ( (TEN8 % var) % 62500L ) / 250;
  633.   USL var_r3 = (TEN8 % var) % 250;
  634.  
  635.   if (sign == PLUS) {
  636.     if (con <= 40) {     /* arctan(1/5)のみ */
  637.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  638.            i = ex.head, j = ex.tail; i <= j; i++) {
  639.         buf      = *(arc+i) + arc_r * TEN8;
  640.         *(arc+i) = buf / co;
  641.         arc_r    = buf % co;
  642.         buf   = ans_r * var_r1;
  643.         quot  = ans_r * var_q;
  644.         quot += 62500L * (buf / va);
  645.         buf2  = 250 * (buf % va) + ans_r * var_r2;
  646.         quot += 250 * (buf2 / va);
  647.         buf   = 250 * (buf2 % va) + ans_r * var_r3 + *(arc+i);
  648.         *(ans+i) += buf / va + quot;
  649.         ans_r     = buf % va;
  650.       }
  651.       for (j = ex.max; i <= j; i++) {
  652.         buf   = ans_r * var_r1;
  653.         quot  = ans_r * var_q;
  654.         quot += 62500L * (buf / va);
  655.         buf2  = 250 * (buf % va) + ans_r * var_r2;
  656.         quot += 250 * (buf2 / va);
  657.         buf   = 250 * (buf2 % va) + ans_r * var_r3;
  658.         *(ans+i) += buf / va + quot;
  659.         ans_r     = buf % va;
  660.       }
  661.     } else {             /* arctan(1/5)以外 */
  662.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  663.            i = ex.head, j = ex.tail; i <= j; i++) {
  664.         buf      = *(arc+i) + arc_r * con_r;
  665.         *(arc+i) = buf / co + arc_r * con_q;
  666.         arc_r    = buf % co;
  667.         buf   = ans_r * var_r1;
  668.         quot  = ans_r * var_q;
  669.         quot += 62500L * (buf / va);
  670.         buf2  = 250 * (buf % va) + ans_r * var_r2;
  671.         quot += 250 * (buf2 / va);
  672.         buf   = 250 * (buf2 % va) + ans_r * var_r3 + *(arc+i);
  673.         *(ans+i) += buf / va + quot;
  674.         ans_r     = buf % va;
  675.       }
  676.       for (j = ex.max; i <= j; i++) {
  677.         buf   = ans_r * var_r1;
  678.         quot  = ans_r * var_q;
  679.         quot += 62500L * (buf / va);
  680.         buf2  = 250 * (buf % va) + ans_r * var_r2;
  681.         quot += 250 * (buf2 / va);
  682.         buf   = 250 * (buf2 % va) + ans_r * var_r3;
  683.         *(ans+i) += buf / va + quot;
  684.         ans_r     = buf % va;
  685.       }
  686.     }
  687.   } else {    /* sign == MINUS */
  688.     if (con <= 40) {     /* arctan(1/5)のみ */
  689.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  690.            i = ex.head, j = ex.tail; i <= j; i++) {
  691.         buf      = *(arc+i) + arc_r * TEN8;
  692.         *(arc+i) = buf / co;
  693.         arc_r    = buf % co;
  694.         buf   = ans_r * var_r1;
  695.         quot  = ans_r * var_q;
  696.         quot += 62500L * (buf / va);
  697.         buf2  = 250 * (buf % va) + ans_r * var_r2;
  698.         quot += 250 * (buf2 / va);
  699.         buf   = 250 * (buf2 % va) + ans_r * var_r3 + *(arc+i);
  700.         *(ans+i) -= buf / va + quot;
  701.         ans_r     = buf % va;
  702.       }
  703.       for (j = ex.max; i <= j; i++) {
  704.         buf   = ans_r * var_r1;
  705.         quot  = ans_r * var_q;
  706.         quot += 62500L * (buf / va);
  707.         buf2  = 250 * (buf % va) + ans_r * var_r2;
  708.         quot += 250 * (buf2 / va);
  709.         buf   = 250 * (buf2 % va) + ans_r * var_r3;
  710.         *(ans+i) -= buf / va + quot;
  711.         ans_r     = buf % va;
  712.       }
  713.     } else {             /* arctan(1/5)以外 */
  714.       for (co = con, va = var, arc_r = 0, ans_r = 0,
  715.            i = ex.head, j = ex.tail; i <= j; i++) {
  716.         buf      = *(arc+i) + arc_r * con_r;
  717.         *(arc+i) = buf / co + arc_r * con_q;
  718.         arc_r    = buf % co;
  719.         buf   = ans_r * var_r1;
  720.         quot  = ans_r * var_q;
  721.         quot += 62500L * (buf / va);
  722.         buf2  = 250 * (buf % va) + ans_r * var_r2;
  723.         quot += 250 * (buf2 / va);
  724.         buf   = 250 * (buf2 % va) + ans_r * var_r3 + *(arc+i);
  725.         *(ans+i) -= buf / va + quot;
  726.         ans_r     = buf % va;
  727.       }
  728.       for (j = ex.max; i <= j; i++) {
  729.         buf   = ans_r * var_r1;
  730.         quot  = ans_r * var_q;
  731.         quot += 62500L * (buf / va);
  732.         buf2  = 250 * (buf % va) + ans_r * var_r2;
  733.         quot += 250 * (buf2 / va);
  734.         buf   = 250 * (buf2 % va) + ans_r * var_r3;
  735.         *(ans+i) -= buf / va + quot;
  736.         ans_r     = buf % va;
  737.       }
  738.     }
  739.   }
  740. }
  741.  
  742. void regular( int mode ) {    /* ans[]を0～TEN8-1の範囲に収める */
  743.   int i;
  744.   long buf, carry;
  745.   long *ans = ex.ans;
  746.   buf   = 0;
  747.   carry = 0;
  748.  
  749.   if (mode == OFF) {    /* 計算中は高速化の為、負の処理をしない */
  750.     for (i = ex.max; i >= 0; i--) {
  751.       buf = *(ans + i) + carry;
  752.       *(ans + i) = buf % TEN8;
  753.       carry      = buf / TEN8;
  754.     }
  755.   } else {              /* 負の処理も行う */
  756.     for (i = ex.max; i >= 0; i--) {
  757.       buf = *(ans + i) + carry;
  758.       *(ans + i) = buf % TEN8;
  759.       carry      = buf / TEN8;
  760.       if ( *(ans + i) < 0) {    /* 負の% 演算の処理 */
  761.         *(ans + i) += TEN8;
  762.           carry--;
  763.       }
  764.     }
  765.   }
  766. }
  767.  
  768. void load( char *fname ) {
  769.   FILE *fp;
  770.   char buf[80];
  771.   int m;
  772.  
  773.   if ( ( fp = fopen(fname, "r") ) == NULL ) {
  774.     fputs( "file not found", stderr );
  775.     usage( !DATA );
  776.   }
  777.   ex.load = ON; /* ロードフラグON */
  778.  
  779.   fgets( buf, 79, fp );  /* まずコメント行をダミーリード */
  780.   fgets( buf, 79, fp );  /* 次に配列に変数群を読み込む */
  781.   if ( sscanf( buf, "%d%ld%ld%ld%d%d%d%d%ld%ld%ld",
  782.                &ex.ver,   /* Version  */
  783.                &ex.sec,   /* 計算時間 */
  784.                &ex.cent,  /* 〃1/100  */
  785.                &ex.keta,  /* 計算桁数 */
  786.                &ex.max,   /* 静的最終 */
  787.                &ex.type,  /* 計算公式 */
  788.                &ex.head,  /* 動的先頭 */
  789.                &ex.tail,  /* 動的最終 */
  790.                &ex.save,  /* 退避flag */
  791.                &ex.con,   /* base^2   */
  792.                &ex.var )  /* arc変数  */ 
  793.       != 11                              /* 変数が足りない   */
  794.       || (ex.max != (int)ex.keta/8 + 2)  /* 関係式がおかしい */
  795.       || ! (ex.save && ex.con && ex.var) /* どれかが0 である */
  796.       ) {                                /* どれでも成立すれば異常 */
  797.       fclose( fp );
  798.       fputs( "Not regular file", stderr );
  799.       exit( 3 );
  800.   }
  801.   get_memory();
  802.  
  803.   for ( m = 0; m <= ex.max; m++ ) {
  804.     fgets( buf, 79, fp );
  805.     if ( sscanf( buf, "%ld%ld", ex.ans+m, ex.arc+m ) != 2 ) {
  806.       fclose( fp );
  807.       fputs( "Not regular file\n", stderr );
  808.       exit( 3 );
  809.     }
  810.   }
  811.   fclose( fp );
  812. }
  813.  
  814. void save( void ) {
  815.   FILE *fp;
  816.   int  m;
  817.   time_t cur_time; /* 現在時刻記録用 */
  818.   char fname[13];  /* 8+1+3+NUL = 13 */
  819.   char buf[80];
  820.  
  821.   time_end();
  822.   sprintf( fname, "%ld.cal", ex.keta );
  823.   if ( ( fp = fopen(fname, "w") ) == NULL ) {
  824.     time( &cur_time );
  825.     fprintf( stderr, "file can't be opened at %s", ctime( &cur_time ) );
  826.     time_set();  /* セーブ失敗でも計算は続行 */
  827.     return;
  828.   }
  829.  
  830.   regular( ON );  /* 既にregular()しているのでこれはズルでない */
  831.   fputs( "PI.C(Ver 2.0)(Normal mode) ここはコメント行（79文字以下）\n", fp );
  832.   sprintf( buf, "%d %ld %ld %ld %d %d %d %d %ld %ld %ld\n",
  833.            ex.ver,   /* Version  */
  834.            ex.sec,   /* 計算時間 */
  835.            ex.cent,  /* 〃1/100  */
  836.            ex.keta,  /* 計算桁数 */
  837.            ex.max,   /* 静的最終 */
  838.            ex.type,  /* 計算公式 */
  839.            ex.head,  /* 動的先頭 */
  840.            ex.tail,  /* 動的最終 */
  841.            ex.save,  /* 退避flag */
  842.            ex.con,   /* base^2   */
  843.            ex.var ); /* arc変数  */
  844.   fputs( buf, fp );
  845.  
  846.   for ( m = 0; m <= ex.max; m++ ) {
  847.     sprintf( buf, "%ld %ld\n", *(ex.ans+m), *(ex.arc+m) );
  848.     fputs( buf, fp );
  849.   }
  850.   fclose( fp );
  851.   time_set();  /* 計算開始時刻を再セット */
  852. }
  853.  
  854. void time_set(void) {  /* 開始時間セット   */
  855.   time( &ex.s_sec );   /* 計算開始時刻(秒) */
  856.   ex.s_cent = clock(); /* 〃1/100 秒       */ 
  857. }
  858.  
  859. void time_end( void ) {  /* 経過時間更新 */
  860.   long    sec, cent;
  861.   time_t  e_sec;
  862.   clock_t e_cent;
  863.  
  864.   time( &e_sec );
  865.   e_cent = clock();
  866.   sec  = (long)difftime( e_sec, ex.s_sec );
  867.   cent = ( e_cent - ex.s_cent + 86400 * (long)CLK_TCK ) 
  868.           % ( 86400 * (long)CLK_TCK );
  869.  
  870.   sec /= 86400;
  871.   sec *= 86400;                  /* 一日未満の秒は一旦切捨て、 */
  872.   sec += (cent / (long)CLK_TCK); /*  cent より求めて再加算する */
  873.   cent *= 100;            /* 1/100 秒単位の数にする為 */
  874.   cent /= (long)CLK_TCK;  /* これでなっているハズ    */
  875.   cent  = cent % 100;     /* 1 秒未満 */
  876.  
  877.   ex.sec  += sec;      /* 計算時間更新 */
  878.   ex.cent += cent;     /* 〃           */
  879.   if (ex.cent >= 100) {
  880.     ex.cent -= 100;
  881.     ex.sec++;
  882.   }
  883.   ex.s_sec  = e_sec;    /* 再開時刻更新(time_setと同じ) */
  884.   ex.s_cent = e_cent;   /* 〃1/100 秒 */ 
  885. }
  886.  
  887. void show_ans( void ) {
  888.   int i = 0;    /* 配列count  */
  889.   int l = 0;    /* 文字列位置 */
  890.   char buf[60];
  891.   char *buf_p = buf;
  892.  
  893.   fprintf( stderr, "計算時間は%ld時間", ex.sec/3600 );
  894.   fprintf( stderr, "%0.2ld分%0.2ld.%0.2ld秒でした\n\n", 
  895.                    (ex.sec/60)%60, ex.sec%60, ex.cent );
  896.   printf("πの小数点以下%d(＋α）桁の値です。\n", ex.keta );
  897.  
  898.   printf( "00000001:%0.1ld.", *ex.ans );
  899.   for ( i = 1; i <= ex.max; i++) {
  900.     sprintf( buf_p, "%0.8ld\0", *(ex.ans + i) );
  901.     buf_p += 8;
  902.     if ( strlen(buf) >= 50 || i == ex.max ) {
  903.       for ( l = 0; l < 50 && buf[l]; l++) {
  904.         putchar( buf[l] );
  905.         if ( l % 10 == 9 )
  906.           putchar( ' ' );
  907.       }
  908.       if (l != 50)
  909.         break;
  910.       printf( "\n%0.8d:  ", (i*8 / 10)*10 + 1 );
  911.       strcpy( buf, buf + 50 );  /* 50桁分詰める */
  912.       buf_p -= 50;              /* 50桁分戻る   */
  913.     }
  914.   }
  915.   printf( "\n計算時間は%ld時間", ex.sec/3600 );
  916.   printf( "%0.2ld分%0.2ld.%0.2ld秒でした\n", 
  917.            (ex.sec/60)%60, ex.sec%60, ex.cent );
  918. }
  919.  
  920. void usage( int mode ) {
  921.   if (mode == DATA) {
  922.     puts("TECNICAL DATA････通常の使い方は  A>RUN386 PI で表示\n");
  923.     puts("TOWNS 2F(初代) 1万桁  10万桁  100万桁");
  924.     puts("ｳｪｲﾄ無調整     96秒   165分   333時間");
  925.     puts("公式別時間     ﾏﾁﾝ    ｶﾞｳｽ    ｼｭﾃﾙﾏｰ  ﾗｻﾞﾌｫｰﾄﾞ");
  926.     puts("               96秒   131秒   141秒   141秒");
  927.     puts("\nセーブファイルについて");
  928.     puts("  サイズ：  　桁数×(1.5～2.1)倍");
  929.     puts("  ファイル名：セーブ時：桁数.CAL（.CALは自動的に付く）");
  930.     puts("              ロード時：拡張子も含め任意の名が可能");
  931.     puts("\n計算再開時にセーブ時間の再変更が可能");
  932.     puts("セーブ指定有りでは、各arctan()終了毎にもセーブする。");
  933.     puts("\n計算進行表示は約１分毎だが途中多少の変動あり");
  934.     puts("\nまた、第２オプションは複数記述が可能");
  935.     exit(0);
  936.   }
  937.   puts("RUN386 PI 100000     ←第１オプション：桁数 or ファイル名");
  938.   puts("          100000.CAL   拡張子必要：ファイルを読み込んで計算再開");
  939.   puts("          (others)     無効オプション： Tecnical data  表示");
  940.   puts("RUN386 PI 100000 G20 ←第２オプション：公式 or セーブ時間間隔(分)");
  941.   puts("                       g:ｶﾞｳｽ m:ﾏﾁﾝ(ﾃﾞﾌｫﾙﾄ) r:ﾗｻﾞﾌｫｰﾄﾞ s:ｼｭﾃﾙﾏｰ");
  942.   puts("                       例はガウス公式で20分おきにセーブする設定");
  943.   puts("約１分毎の計算進行表示時、 Ctrl-C による中断が可能である");
  944.   exit(0);
  945. }
  946.